#list the libraries needed
necessary=c("adehabitat","SDMTools","rgdal","sp")

#check if library is installed
installed = necessary %in% installed.packages()

#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

### This script will reformat the BVG ASCII to 1:2 million

in.dir = '/home1/99/jc152199/TopoGapAnalysis/BVG/'
out.dir = '/home1/99/jc152199/TopoGapAnalysis/ASCIIFiles/'

### Read in template ASCII and relate information

base.asc = read.asc(paste(in.dir,'bvg.asc',sep=""))
relates = read.csv(paste(in.dir,'bvg.csv',sep=""),header=T)

# Create an expanded data frame of row/column positions for the ASCII files

base.pos = as.data.frame(which(is.finite(base.asc), arr.ind = T))

base.pos$east = getXYcoords(base.asc)$x[base.pos$row]
base.pos$north = getXYcoords(base.asc)$y[base.pos$col]

pnts = data.frame(east=base.pos$east, north=base.pos$north)

# Read in other topo ASCII files

slope.asc = read.asc(paste(out.dir,'slope80.asc',sep=""))
aspect.asc = read.asc(paste(out.dir,'aspect80.asc',sep=""))
coastdist.asc = read.asc(paste(out.dir,'coast_dist80.asc',sep=""))
dem.asc = read.asc(paste(out.dir,'dem80.asc',sep=""))

# Extract data from the above four ASCII's based on the east/north combos from base.pos

slope.values = extract.data(pnts,slope.asc)
aspect.values = extract.data(pnts,aspect.asc)
coastdist.values = extract.data(pnts,coastdist.asc)
dem.values = extract.data(pnts,dem.asc)

# Add the value from the base.asc 

base.pos$value = base.asc[cbind(base.pos$row,base.pos$col)]

# Still need to change base.pos$BVG to the 1:2 million level

data1 = merge(base.pos,relates,by='value')#,all.x=T)

# Create the 5 new re-sized ASCII files & print images

slope.asc = base.asc; slope.asc[,] = NA
slope.asc[cbind(base.pos$row,base.pos$col)]= slope.values
png('slope.png')
  image(slope.asc)
dev.off()

aspect.asc = base.asc; aspect.asc[,] = NA
aspect.asc[cbind(base.pos$row,base.pos$col)]= aspect.values
png('aspect.png')
  image(aspect.asc)
dev.off()

coastdist.asc = base.asc; coastdist.asc[,] = NA
coastdist.asc[cbind(base.pos$row,base.pos$col)]= coastdist.values
png('coastdist.png')
  image(coastdist.asc)
dev.off()

dem.asc = base.asc; dem.asc[,] = NA
dem.asc[cbind(base.pos$row,base.pos$col)]= dem.values
png('dem.png')
  image(dem.asc)
dev.off()

BVG.asc = base.asc; BVG.asc[,] = NA
BVG.asc[cbind(data1$row,data1$col)]= data1$BVG_2M
png('BVG.png')
  image(BVG.asc)
dev.off()

# Backup the new ASCII files

write.asc(slope.asc,paste(out.dir,"slope80WTplusbuffer.asc",sep=""))
write.asc(aspect.asc,paste(out.dir,"aspect80WTplusbuffer.asc",sep=""))
write.asc(coastdist.asc,paste(out.dir,"coastdist80WTplusbuffer.asc",sep=""))
write.asc(dem.asc,paste(out.dir,"dem80WTplusbuffer.asc",sep=""))
write.asc(BVG.asc,paste(out.dir,"BVG2M80WTplusbuffer.asc",sep=""))

# Read in topo ASCII files

slope.asc=read.asc(paste(out.dir,"slope80WTplusbuffer.asc",sep=""))
aspect.asc=read.asc(paste(out.dir,"aspect80WTplusbuffer.asc",sep=""))
coastdist.asc=read.asc(paste(out.dir,"coastdist80WTplusbuffer.asc",sep=""))
dem.asc=read.asc(paste(out.dir,"dem80WTplusbuffer.asc",sep=""))
BVG.asc=read.asc(paste(out.dir,"BVG2M80WTplusbuffer.asc",sep=""))

# Now create a subset of base.pos for the locations of all reptile samples

surveylocs = read.csv(paste('/home1/99/jc152199/TopoGapAnalysis/cjs_Rsurveylocs.csv',sep=""),header=T)

survey.pnts = cbind(surveylocs$east,surveylocs$north)

# Use extract.data command to get info from the ASCII files for all survey locations

survey.slope = extract.data(survey.pnts, slope.asc)
survey.coastdist = extract.data(survey.pnts, coastdist.asc)
survey.aspect = extract.data(survey.pnts, aspect.asc)
survey.dem = extract.data(survey.pnts, dem.asc)
survey.BVG = extract.data(survey.pnts, BVG.asc)
survey.BVGRFonly = extract.data(survey.pnts, BVGRFonly.asc)

# Now trying to create an ASCII file of rainforest patch size

BVGbinary.asc = BVG.asc
BVGbinary.asc[which(BVG.asc>=7)]=1
BVGbinary.asc[which(BVG.asc<=8)]=0
BVGRFonly.asc = BVG.asc
BVGRFonly.asc[which(BVG.asc>7)]=0

# Run ConnCompLabel on Binary Matrix to assign unique numbers to individual patches

RFpatches.asc = ConnCompLabel(BVGbinary.asc)

# Produce histograms of whole region

png(filename=paste(out.dir,"BVGRFonlynewest.png",sep=""), units="cm", res=100, height=40, width=80, bg="white")

  #create the plot

  histbackback(slope.asc, survey.slope, probability=TRUE, brks=seq(0,70,by=2), xlim=c(-.2,.2))
  histbackback(aspect.asc, survey.aspect, probability=TRUE, brks=seq(-1,374,by=15))
  histbackback(coastdist.asc, survey.coastdist, probability=TRUE, brks=seq(0,90000,by=5000), xlim=c(-.0004,.0004))
  histbackback(dem.asc, survey.dem, probability=TRUE, brks=seq(0,1600,by=100))
  histbackback(BVG.asc, survey.BVG, probability=TRUE, brks=seq(0,100,by=1), xlim=c(-.4,.4))
  histbackback(BVGRFonly.asc, survey.BVGRFonly, probability=TRUE, brks=seq(-.5,8.5,by=1), xlim=c(-1,1))
  
  #Turn off device driver (to flush output)

  dev.off()
